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1.  Introduction 

It  is  becoming  increasingly  important  to  use  computer-aided  device  simulators  in  the 
design  and  fabrication  processes  of  semiconductor  devices  designing  as  semiconductor 
devices  continue  to  decrease  in  size  toward  the  deep  submicron  regime.  The  development 
of  devices  involves  several  iterations  of  trial  and  error  in  fabrication  imtil  a  specified  goal 
in  terms  of  design  conditions  is  reached.  The  application  of  device  modeling  can  provide 
an  inexpensive  way  to  analyze  and  design  the  semiconductor  devices  before  expensive 
device  processing.  Since  traditional  equivalent  circuit  models  and  close-form  analytical 
models  cannot  always  provide  consistently  accurate  results  for  all  modes  of  operation  of 
today’s  small  devices.  This  has  meant  that  there  has  been  a  greater  demand  for  models 
capable  of  increasing  our  understanding  of  how  these  devices  operate  and  capable  of 
predicting  accurate  quantitative  results. 

To  simulate  a  device,  we  solve  a  transport  equation  coupled  with  Poisson  equation 
The  accuracy  of  a  simulation  is  usually  determined  by  how  accurately  carrier  transport  is 
described.  Generally,  the  more  sophisticated  the  approach,  the  heavier  the  computational 
burden,  so  it  is  important  to  choose  an  adequate  approach  for  the  device  under  study.  In 
the  past,  the  study  of  electric  behavior  in  a  semiconductor  device  has  been  based  on  the 
drift-diffiision  equation.  The  drift-diffusion  model  is  a  low-order  approximation  of  the 
Boltzmann  transport  equation,  it  implies  that  mobility  of  the  carrier  is  only  a  ftmction  of 
the  local  electrical  field  and  it  does  not  take  account  of  the  non-stationary  characteristics 
such  as  carrier  heating  and  velocity  overshoot  [1,2].  The  application  of  this  model  is 
limited  to  devices  where  the  spatial  variation  of  the  electric  field  is  not  very  large. 
However,  in  modem  devices,  whose  size  is  in  the  deep  submicron  region,  the  non- 
stationary  phenomena  are  becoming  more  important.  As  a  result,  the  drift-diffusion 
model  is  no  longer  applicable  [3-6]. 

Monte  Carlo  simulation  has  been  widely  used  for  analyzing  carrier  transport  in  bulk 
semiconductors  [7,8,9].  This  method  tracks  the  momentum  and  position  of  an  ensemble 
of  carriers  as  they  move  through  a  device  imder  the  influence  of  an  electric  field  and 
random  scattering  forces.  Random  numbers  are  chosen  to  determine  the  time  between 
collisions,  the  type  of  scattering  events  encountered,  and  the  direction  of  the  carrier  after 
scattering.  Monte  Carlo  solution  of  the  Boltzmann  transport  Equation  can  provide  a  more 
detailed  description  of  carrier  transport.  This  is  because  band  structure  and  various 
scattering  mechanisms  can  be  taken  into  consideration  [10,1 1].  However,  the  noise  in  the 
solution  and  the  enormous  amounts  of  computation  time  required  make  this  model 
impractical  for  device  design. 

Another  device  simulation  method  is  based  on  the  hydrodynamic  model,  which  is 
obtained  by  taking  the  first  three  moments  of  the  Boltzmann  Transport  Equation  [12-15]. 
These  are  the  carrier  continuity  equation,  the  momentum  balance  equation,  and  the 
energy  balance  equation.  Instead  of  solving  for  the  carrier  distribution  function,  the 
quantities  of  interest  are  calculated  directly  through  the  hydrodynamic  equations.  The 
hydrodynamic  model  has  many  advantages  over  ftie  drift-diffusion  model,  such  as  its 
ability  to  treat  high-field,  nonstationary,  and  hot-electron  effects.  It  can  also  model  multi¬ 
dimensional  devices  without  the  excessive  computational  cost  of  the  Monte  Carlo 
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simulations.  But  it  has  been  known  to  give  unphysical  solutions  in  some  cases,  e.g.,  the 

spurious  velocity  overshoot  spike  in  ‘ballistic’  diode  [16].  The  moments 

equations  by  themselves  do  not  form  a  closed  set  of  equations.  The  momentum  and 
energy  relaxation  times  are  needed  and  can  only  be  imported  from  outside  of  the  system. 
These  relaxation  times  are  supplied  from  Monte  Carlo  simulations  or  simply  assumed  to 
be  constant  [17,18]. 

Recently,  a  hydrodynamic-balance-equations  model  of  carrier  transport  in 
semiconductor  devices  based  on  the  Lei-Ting  balance  equations  has  been  developed 
[19,20,21].  The  Lei-Ting  balance  equations  have  been  successfully  applied  to  many  types 
of  semiconductor  microstructures,  including  studies  of  nonequilibrium  phonon, 
nonstationary  and  high  frequency  transport  [22,23].  Unlike  other  hydrodynamic  equation 
based  approaches  to  device  modeling,  where  the  various  relaxation  rates  are  imported 
from  Monte  Carlo  calculations  or  simply  assumed  to  be  constant  [24,25,26],  the  Lei-Ting 
hydrodynamic  balance  equations  approach  includes  scattering  in  the  form  of  frictional 
force  functions  due  to  electron-impurity  and  electron-phonon  interaction  and  an  energy 
loss  function  due  to  electron-phonon  interaction.  These  quantities  are  calculated  within 
the  simulation  process  itself,  as  functions  of  the  electron  drift  velocity,  electron 
temperature,  as  well  as  the  electron  density,  without  an  outside,  separate  Monte  Carlo 
procedure  [21,27,28].  Thus,  besides  the  usual  advantages  of  traditional  hydrodynamic 
simulation  approaches,  the  present  method  enjoys  the  added  convenience  of  self- 
contained  treatment  of  scattering. 

This  document  reports  on  a  systematic  implementation  of  the  Lei-Ting  hydrodynamic 
balance  equations  as  a  sophisticated,  versatile  device  simulation  package,  capable  of  ID, 
2D  device  modeling  tasks  encountered  by  device  designers  today.  In  addition  to  steady- 
state  modeling,  transient  device  simulations  based  on  the  new  hydrodynamic  model  are 
also  described.  Without  any  complicated  mathematics,  a  decoupled  method  with  a 
relatively  large  time  step  has  been  applied  to  the  transient  simulation.  The  time 
discretization  algorithm  for  our  transient  simulator  is  based  on  the  Crank-Nicolson 
method  for  time  discretization,  and  for  the  algorithm  of  the  time  step  selection  we  uses 
the  local  error  to  determine  the  size  of  each  time  step. 

This  report  is  organized  as  follows.  In  chapter  2,  we  give  a  brief  review  of  the 
commonly  used  models  of  semiconductor  transport,  such  as  Monte  Carlo  method,  drift- 
diffusion  model  and  hydrodynamic  model.  A  nonparabolic  multivalley  Lei-Ting  balance 
equation  is  given  in  Chapter  3.  Hydrodynamic  balance  equations  for  a  single  valley 
semiconductor  have  been  described  in  Chapter  4.  In  Chapter  5,  the  discretization  of  the 
hydrodynamic  balance  equations  is  performed.  The  Box  Integration  Method  is  used  for 
spatial  discretization  and  the  Crank-Nicolson  implicit  scheme  is  used  for  time 
discretization.  In  Chapter  6,  we  present  some  simulation  results  of  the  application  of  our 
hydrodynamic  model.  Discussion  and  conclusion  are  given  in  Chapter  7. 
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2.  Models  of  Semiconductor  Transport 

2.1  Boltzmann  Transport  Equation  and  Monte  Carlo  Method 


To  simulate  carrier  motion  in  a  submicron  semiconductor  device,  one  has  to  solve  the 
Boltzmann  transport  equation  and  the  Poisson  equation  simultaneously.  The  Boltzmann 
transport  equation  is  given  by: 


dt 


rJ  ^  kJ 


(2.1) 


where  /  =  f(f,k,t)  is  the  distribution  function,  t  is  the  time  variable,  v  is  the  carriers 

velocity,  F  is  the  force  on  the  carriers.  The  collision  term  on  the  right-hand  side  of 
equation  (1)  includes  all  the  scattering  mechanisms,  and  it  is  given  by 

(f)^  (2-2) 

Here,  W{k,k')  is  the  transition  probability,  and  W{k,k')dV^dt  is  the  conditional 

probability  of  the  transition  from  the  state  k'  in  dV^.  in  time  dt  given  the  an  electron  is 

initially  in  state  k  and  the  state  k'  is  empty.  A  direct  solution  to  the  Boltzmann  transport 
equation  coupled  to  a  self-consistent  electric  field  pattern  for  any  realistic  semiconductor 
devices  is  rattier  difficult  because  the  equation  to  be  solved  is  a  very  complicated  integro- 
differential  equation  with  seven  independent  variables. 

The  most  popular  technique  for  solving  the  Boltzmann  transport  equation  is  the 
Monte  Carlo  method.  This  method  tracks  the  position  and  momentum  of  an  ensemble  of 
particles  as  they  move  through  a  device  under  the  influence  of  the  electric  field  and 
random  scattering  forces.  Random  numbers  are  chosen  to  determine  the  time  between 
collisions,  the  type  of  scattering  events  encountered,  and  the  direction  of  the  carrier  after 
scattering.  The  processes  repeated,  typically  between  10"*  and  10®  times,  to  simulate  the 
carrier  path  through  the  device.  For  the  time  dependent  problems,  the  sample  ensemble 
must  be  sufficiently  large  to  accurately  represent  the  actual  electron  gas. 

The  Monte  Carlo  simulation  can  provides  more  accurate  simulations  of  carrier 
transport  in  a  semiconductor  device,  it  does  have  some  limitations  that  should  be 
mentioned.  First,  the  intensive  computer  time  required  and  the  statistical  noise  associated 
with  it  limit  this  method’s  application.  Although  the  statistical  noise  decreases  as  the 

number  of  simulated  carriers,  Nsim,  increases,  but  it  only  drops  as  ,  so 

unreasonably  large  number  of  carriers  would  have  to  be  simulated.  Second,  Monte  Carlo 
method  is  not  well  suited  for  dealing  with  the  low-field  region,  barrier  region  and  carrier 
generation-recombination  processes. 
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2.2  The  Drift-Diffusion  model 


Currently,  most  device  modeling  packages  are  based  on  the  drift-diffusion  (DD)  model, 
which  consists  of  the  basic  semiconductor  equations 


where 


ot  e 

V.(eV(|>)  =  e(i\r^-n), 


(2.3) 

(2.4) 


J„=eD„Vn-en\i„V^.  (2.5) 

Equation  (2.3)  is  the  carrier-continuity  equation,  and  (2.4)  is  the  Poisson  equation.  Here, 
is  the  electron  current  density,  n  is  the  electron  density,  is  the  electron  mobility,  U 
is  the  net  recombination  rate  per  unit  volume,  ^  is  ftie  electrical  potential,  andiVo  is  the 
donor  density.  Diffusivity  and  mobility  are  assumed  to  be  related  by  the  Einstein  relation 


(2.6) 


where  T  is  the  lattice  temperature.  The  mobility  is  related  to  the  momentum  relaxation 
time  and  effective  mass  as  follows: 


ez. 


(2.7) 


Here  and  are  momentum  relaxation  time  and  effective  mass  of  electron. 

The  drift-diffusion  model  can  be  shown  to  arise  from  the  second-moment 
approximation  to  the  Boltzmann  transport  equation.  Since  the  drift-diffusion  model  are 
based  on  the  assumptions  that  carriers  are  in  local  equilibrium  with  the  lattice,  and  the 
mobility  is  a  function  of  the  local  electric  field,  the  applicability  of  this  model  is  limited 
to  devices  where  the  spatial  variation  of  the  electric  field  is  not  very  large.  Due  to  rapid 
variations  of  the  electric  field  in  today’s  device,  these  assumptions  are  not  always  valid, 
and  hence  the  transport  coefficients  cannot  be  considered  to  be  dependent  only  on  the 
local  electrical  field.  The  drift-diffusion  model  fails  to  predict  phenomena  such  as  carrier 
heating  and  velocity  overshoot.  Velocity  overshoot  can  be  explained  if  we  recognize  that 
the  mobility  is  more  closely  related  to  the  local  carrier  energy  than  the  local  electric  filed. 
Velocity  overshoot  occurs  because  at  the  beginning  of  the  high  field  region  the  mobility 
is  at  its  low  field  value.  As  the  carrier  energy  rises,  the  mobility  drops  until  the  saturated 
velocity  is  reached.  Velocity  overshoot  is  a  non-local  effect  and  also  is  referred  to  as  an 
off-equilibrium  effect,  because  carriers  do  not  have  an  opportunity  to  establish 
equilibrium  with  the  local  electric  field  in  submicron  devices.  Velocity  overshoot  can 
improve  the  performance  of  small  device  because  the  average  carrier  velocity  can  exceed 
the  steady  state  limit. 


2.3  The  Hydrodynamic  model 


The  hydrodynamic  model  consists  of  a  set  of  equations  expressing  the  conservation  of 
charge,  momentum  and  energy  for  each  species  of  carriers  [4].  These  equations  are 
derived  by  taking  the  first  three  moments  of  the  Boltzmann  transport  equation.  The 
hydrodynamic  equations  for  electrons  can  be  written  as  follows: 
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(2.8) 

^  ( J  \  Jr 

'  =cD,Vn  +  enn,(-5-Vj;-V(|)), 
atynj  e 

(2.9) 

5(«w„)  ^  s  -i  =  dw 

(2.10) 

1 

11 

> 

CO 

> 

(2.11) 

where  Tgis  the  electron  temperature,  is  the  electron  momentum  relaxation  time,  w„  is 

the  electron  mean  energy,  E  =  -V(|)  is  the  electric  field.  Here,  the  diffusion  coefficient  is 
given  by  the  generalized  Einstein  relation 

(2.12) 

e 

where  is  the  electron  temperature. 


The  electron  energy  flow  S„  is  given  by 

=  -K„vr^  -(w„  +kgT;) ,  (2.13) 

The  thermal  conductivity  k„  is  related  to  the  mobility  by  the  Wiedmann-Franz  law: 

K.  =(i+c)env.X^fT„  (2.14) 

2  e 

where  c  is  an  adjustable  parameter,  usually  taken  to  be  -1. 

The  electron  mean  energy  w„  is  defined  as 

,  (2.1 5) 

m„  the  electron  effective  mass,  and  v„  the  mean  velocity  of  electron. 

The  collision  terms  are  modeled  by  the  relaxation-time  approximation.  The  collision 
term  in  (2.10)  is  expressed  as: 


dw  W  —  Wn 

n  \  _  n  0 

V  ~ 


(2.16) 


where  is  the  energy  relaxation  time  for  the  electron. 

Two  parameters  are  needed  in  applying  the  hydrodynamic  model:  the  momentum 
relaxation  time  and  the  energy  relaxation  time.  A  common  approach  is  to  extract  these 
parameters  jfrom  bulk  homogeneous  Monte  Carlo  simulations.  An  alternative  method  is  to 
use  various  empirical  mobility  models,  such  as  [29]: 


^  ^^0  .  y=_J£o_^ 

1  +  Y(w-Wo)’  ’ 


(2.17) 


where  is  the  electron  saturation  velocity,  po  is  the  low-field  mobility,  and  is  the 
low-filed  energy  relaxation  time  in  the  homogeneous  case. 

Although  hydrodynamic  models  include  more  physics  than  does  the  drift-diffusion 
model,  serious  errors  may  occur  in  practice,  e.g.,  the  spurious  velocity  overshoot  spike  in 


ti^  -n-ri^  ‘ballistic’  diode.  This  spurious  peak  can  be  eliminated  by  varying  the  value 
of  the  thermal  conductivity  in  the  Wiedmann-Franz  law.  However,  the  value  of  the 
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thermal  conductivity  depends  on  the  doping  distribution  and  on  the  applied  voltages, 
therefore  it  has  to  be  adjusted  to  fit  the  individual  devices. 
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3.  The  Lei-Ting  Hydrodynamic  Balance  Equations  for  a  Nonparabolic  Multivalley 
Semiconductor 


The  Lei-Ting  balance  equation  approach  for  hot  electron  transport  in  semiconductors 
is  based  on  the  Heisenberg  equations  of  motion  for  the  total  physical  momentum,  the 
total  energy  and  the  population  of  carriers  in  each  energy  valley,  and  the  statistical 
average  v^ith  respect  to  an  initial  density  matrix  having  a  lattice  wave-vector  shift,  an 
electron  temperature,  and  a  chemical  potential  for  each  energy  valley  as  parameters. 

Here,  we  consider  a  semiconductor  system  with  ann-valley  band  structure  in  a  electric 

field  E  and  the  energy  band  &^{k)  has  a  general  spectrum  functions  in  the  £  space.  The 

equations  for  momentum,  energy  and  carrier  numbers  balance  in  the  ath  valley  are 
written  in  the  form  [30]: 

‘':Wv,)  =  Ar.e£-K,+A”+A; 


dt 


d 


+  Z  A +  E  A 


ab 


b(^a) 


b(^a) 


+  Z  A 

b(^a) 


ab 


^(JVA)  =  N,eE-v, -W^-YW‘;-YK, 


dt 

b{i^a)  b{i^a) 


b{^a) 


bi^a) 


(3.1) 


(3.2) 


(3.3) 


Here,  Na  is  the  average  number,  is  the  average  velocity,  is  the  average  energy, 
is  the  ensemble-averaged  inverse  effective  mass  tensor.  A"  is  the  frictional  acceleration 
contributed  from  the  intravalley  scattering,  A“*is  the  frictional  acceleration  contributed 
from  the  intervalley  scattering,  A““  is  the  frictional  acceleration  contributed  from  the 
electron  intravalley-intravalley  Coulomb  scattering,  is  the  energy  loss  rate  due  to 

the  intravalley  electron-phonon  scattering,  W^is  the  energy  loss  rate  due  to  the 

intervalley  scattering,  the  energy  loss  rate  due  to  the  electron  mtravalley- 

intravalley  Coulomb  scattering,  and  X"*  is  the  rate  of  numbers  change  of  the  carriers  in 
the  ath  energy  valley.  The  statistic  averages  of  the  relevant  quantities  are  defined  by  the 
following 


K 


= (^)yT(8,  (^)  -  p  J  /  rj , 

k 

a  k 

= (^)]/[(^.  (^)  -  p  j  /?;  ] . 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


(3.8) 
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where  /(x)  =  l/[exp(A:)  +  l]  is  the  Fermi-Dirac  function,  is  the  average  lattice 
momentum,  is  the  chemical  potential,  is  the  carrier  temperature,  and  the  velocity 
function  of  carriers  is  given  by 

=  (3.9) 

The  average  drift  velocity  of  the  whole  system  is  given  by 

=  (3.10) 

a 

where  N  is  the  total  carrier  number  in  this  system. 

The  intravalley  frictional  acceleration  due  to  the  electron-impurity  scattering  is  given 
by 

K  =  2]  \u{q)f  [v„  {k)  -v^{k  +  ^)]8k  (^)  -  (^  +  ^)] 

=  I 


/(e,(^),7;)-/(6,(^+g),rj 


(3.11) 


and  the  intravalley  frictional  acceleration  due  to  the  electron-phonon  scattering  is  given 
by 

Kp=^  -45^2 1^ (k)-v^(k+  ^)]6[s^ (^)  - 8^ (^  +  q)\ 


(3.13) 


^  /(s,(nrj-/(e,(£  +  g).r,)  r _  r  e,(^)-8,(^  +  g)^  _  (3  J2) 

[9.s„(^)-£„(^  +  ^)]  ^  _V^bT  )  t  kgT^  J_ 

The  energy  loss  rate  contributed  from  the  intravalley  scattering  is  expressed  as 
=  47c2  ^)r  |gaa  ik,  q)^  Q^,;,6[8„  (^)  -  8„  (^  +  ^)  -  ] 

k,qX 

X  fi^a(k)Ja)-f{^a(k  +  q),Ta)\(  _  (  S,  (^)  -  8,  (^  +  g)  ^1  (3^13) 

The  intervalley  frictional  acceleration  due  to  the  electron-impurity  scattering  A“f  is 
given  by 

A"f  =  -27mi  J]|u(g)f  {k, qf  v„ (^)5[8^ {k)  -  s* (^  +  q)] 

k,q 

X  [/(8.  {hTa  )  -  /(8,  {k  +  q),T,)-\ .  (3.14) 

And  the  intervalley  frictional  acceleration  due  to  electron-phonon  scattering  A°*  is  given 
by 


AJ  =-AnY,\M{q,xf  g^{k,q)  v„(^)At(^,g,n^J. 


(3.15) 


The  frictional  acceleration  contributed  from  the  electron  intravalley  Coulomb  scattering 
A“*  is  given  by 

Af  =  -2X  K  ^^f\saa  (^.  4i\\sbb  {k\qi\[v^  (k)  -vjk  +  q)]Y‘'^  (k,k\q)  ■  (3.16) 

k,k\q 

The  energy-loss  rate  due  to  the  intervalley  electron-phonon  scattering  is  expressed  as 
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(3.17) 

k^q^X 

and  the  intravalley  Coulomb  scattering  energy-loss  rate  is  given  by 

Ke  =  -2  -8„(^  +  ^)]Y“*(^J',9)  (3.18) 


k,k  yq 


The  rate  of  numbers  change  of  the  carriers  in  the  ath  energy  valley  due  to  the  electron- 
impurity  scattering  is  given  by 

X“f  =  -27t«, 2  \u{qf  (k,q)^ 6[8,  (^)  -  8^ {k  +  q)\f(&,  (k),  )  - /(Sj  (k  +  q),  Tj )] ,  (3. 1 9) 

k>q 

and  the  rate  of  numbers  change  of  carriers  due  to  the  electron-phonon  scattering  is  given 
by 

K  =  At(r,?,n„).  (3.20) 

k,g,X 

Here, 

+|/(8,(^),r,)-/(8,(^+9),r,) 

and 

Y°\k,k\d)  =  27r5[8„  (^) -E^{k +q)  +  B^{k'  +q)- s  ^  (£')] 

X  [9.s„(^)-s^(^  +  9)]  ef'"  [9,6j(^')-Sj (*’  +  ?)] 

x[/(£„(^),rj-/(s„(^+a7;)][/(s,(a?;)-/(8,(^  +  a7;)] 


n 

-n 

\(^)  ej(^  +  9)'| 

<  ) 

V  J 


-n 


U^jk  +  q)  ejk)'^ 

\  kgTi,  kgT^  j 


5[s,(^)-8,(*  +  9)  +  fin^],  (3.21) 


(  s^ik)-s„ik  +  q)^ 

-  n 

(  Bgik')- Zgiic'  +  ^)"| 

1  kgT^  ) 

(3.22) 


and  the  Bose  distribution  function  is  expressed  as 

nix)  =  ^ 


(3.23) 


exp(x)  - 1 

In  the  above  equations,  «i  is  the  impurity  density,  is  the  phonon  frequency  of  the  A,th 
branch  with  wave  vector  q,  the  intravalley,  uiq)  is  the  electron-impurity  potential, 
Miq,X)  is  the  electron-phonon  matrix  element,  v^(^)  is  the  electron-electron  Coulomb 

potential  in  the  plane-wave  representation,  and  gaaik,q)  andg„j(^,^)  are,  respectively, 
the  intravalley  and  intervalley  form  factors  related  to  the  wave  functions  of  the  oth  and 
hth  energy  valley, 


with 


= 


1  +  azik) 

1/2 

asik) 

\  +  2aj&ik)_ 

5 

l  +  2asik)_ 

(3.24) 

(3.25) 
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For  intravalley  scattering,  we  consider  the  intravalley  carrier  screening  in  the  random- 

phase  approximation  (RPA),  ef  ^  (^,©)  is  the  RPA  dielectric  function  of  the  ath  valley. 

For  a  Kane-type  nonparabolic  energy  band,  the  relationship  between  the  energy  and 
the  wave  vector  can  be  written  as 


s,(k)=^ 

2a 


(3.26) 


where  a  is  the  nonparabolic  parameter.  The  parabolic  energy  dispersion  of  the  ath  valley 
is  given  by 


,(P) 


(k)  = 


iK-Kf  ■  ,  iK-Kf 


2m„ 


2m, 


ay 


2m„. 


(3.27) 


where  Umax,  \lmay,  and  Vmaz  are  theXjy,  andz  components  of  the  inverse  effective  mass 


tensor,  respectively,  and  indicates  the  position  in  the  Brillouin  zone  of  the  center  of 
the  ath  valley. 

The  scattering  mechanisms  which  have  been  considered  may  include  intravalley 
electron-impurity  scattering,  intravalley  electron-phonon  scattering  (acoustic  and  optic), 
intervalley  electron-phonon  scattering  (acoustic  and  optic),  and  electron-electron 
interaction. 

For  Si,  the  electrons  which  contribute  to  transport  are  those  in  the  six  equivalent 
valleys  which  lie  along  the  <100>  directions,  as  shown  in  Fig  3.1. 


Fig.  3.1  Constant  energy  surfaces  for  the  conduction  band  of  silicon. 

There  are  two  types  of  intervalley  scattering  in  Si:  g-type  (between  parallel  valleys) 
and  y-type  (between  perpendicular  valleys).  So  that,  the  squared  matrix  elements  forg- 
type  and y-type  scatterings  are  represented  by 
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I  .2 


and 


2Cl^^d 


I  l2 


(for  a=l,  6=2), 


(for  a=l,  6=3, 4, 5, 6). 


For  the  acoustic  intravalley  scattering,  the  squared  matrix  elements  is  given  by 


Kr= 


nE^q 


2v,d 


(3.28) 

(3.29) 


(3.30) 
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Table  1.  Physical  parameters  for  Silicon[30]. 


Mass  density 

d 

2.329  g/cm" 

Longitudinal  soimd  velocity 

Vs 

9.04x10^  cm/s 

Longitudinal  effective  mass 

Ml 

0.916  mo 

Transverse  effective  mass 

rrit 

0.19  mo 

Dielectric  constant 

K 

11.7 

Longitudinal  optical  phonon 

h(OQ 

0.063  eV 

energy 

Acoustic  deformation  potential 

El 

9.2  eV 

Nonparabolicity  parameter 

a 

0.5  eV‘ 

Intervalley  scattering 


Equivalent  temperature  (K)  Coupling  constant  (xlO®  eV/cm) 


/-type 

210 

0.15 

500 

3.4 

600 

4.0 

g-type 

D, 

A 

140 

0.5 

210 

0.8 

700 

3.0 

4.  Hydrodynamic  Balance  Equations  for  a  Single  Valley  Semiconductor 


4.1  Lei-Ting  Balance  Equation  for  a  Single  Parabolic  Energy  Band 


The  Lei-Ting  balance  equations  for  a  single  parabolic  band  system  can  be  summarized  as 
follows  [19,21]: 


^.iv.(J)  =  -cr, 

at  e 

(4.1) 

dv  2  V«  e  p  / 

dt  3  mn  m  mn 

(4.2) 

^  +  v.Vm  =  --mV-v-F-v-/-V-Q 

8t  3  J 

(4.3) 

V^<t)  =  --(V  -«), 

8  ^ 

(4.4) 

where  n  is  the  statistical  average  of  the  electron  number  density,  v  is  the  average 
velocity  of  the  electrons,  U  is  the  net  recombination  rate  per  imit  volume,  u  is  the 

average  kinetic  energy  of  the  relative  electrons.  The  symbol  /  denotes  the  frictional 
force  experienced  by  the  electrons  due  to  impurity  and  phonon  scattering,  W  is  the 
energy-loss  rate  of  the  electron  system  to  the  phonon  system,  is  the  heat  flow.  The 
symbols  e  and  m  are  the  electron  charge  and  effective  mass,  s  is  the  static  background 
dielectric  constant,  is  the  net  doping  concentration. 

The  average  local  kinetic  energy  density  of  the  relative  electrons  is  given  by 

«(^)  =  2yej/„[(e.--ji(«))/i.2;(R)]  .  (4.5) 

k 

The  local  chemical  potential  \i(R)  is  related  to  the  local  electron  density  n(R)  via  the 
relation 

n(*)  =  2y/o[(e,-n(^))/J:,r.(5)]  ,  (4.6) 

k 

where  is  the  band  energy  of  the  electron,  and  fo  is  the  Fermi-Dirac  function. 

With  the  help  of  the  Fermi  integrals,  the  local  carrier  number  and  the  internal  energy 
densities  can  be  expressed  as 

«  =  \K‘r,N,(J'.V  y&,). 


where 


(4.7) 

(4.8) 
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and  the  effective  density  of  states  in  the  conduction  band  is 


2rdi 


The  Fermi  integral  of  order  j  is  given  by 


(4.10) 


r(y+l)'’o  exp(y-x)  +  l 

At  room  temperature,  the  variable  jc  in  the  Fermi  integral  (4.1 1)  is  much  less  than  zero, 
so  that  die  Fermi  integral  can  be  approximated  by  its  nondegenerate  form 

^(x)  =  exp(x). 


for  not  too  high  densities. 
The  resistive  force  is 


/  =  «/! 

9 


Uiq) 


+2z|m(^x)|  ^n2(^,«)o+^^)x 

qX 


n  ( — —)-n  ( - - — ) 


and  the  energy  loss  rate  is 

W=2X\M(^fn^  n2(?,(»o  +  )  X  ] (^)  - ^ ) 

B 


(4.12) 


(4.13) 


(4.14) 


where  ©o  =q-v{R),  %(x)=(e^-l)  ^  is  the  Bose-Einstein  factor,  «/  is  impurity  density, 
is  the  phonon  frequency  of  wave  vector  q  and  branch  index  X,  17(^)  is  the  Fourier 
transform  of  the  electron-impvirity  interaction  potential,  M{q,X)  is  the  electron-phonon 
coupling  matrix  element,  112  (^,®)  is  the  imaginary  part  of  the  density-density 
correlation  function  of  electrons  that  can  be  obtained  within  the  random-phase 
approximation  (RPA)  or  beyond.  Note  that  /  and  W  depend  on  the  position  vector  R 
through  the  quantities  w(.R) ,  Tg{R),  and  v(^) . 


4.2  Scatting  Mechanisms 

The  total  resistive  force  density /is  composed  of  two  parts,  due  to  impurity  scattering  and 
phonon  scattering: 

f  =  A.p+fpH  .  (4-15) 

which  in  the  case  of  a  three-dimensional  bulk  semiconductor  are  given  as 

u = ,  (4.16) 

=-^^f_dc[dq'Z\qM(qXfxP,(_q.n,^ 
where  =  ^vx ,and  (^^  =  1/ kgTg)  , 


(4.17) 
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=  +®))-  (4.18) 

The  energy  loss  rate  is  determined  by  inelastic  electron-phonon  scattering 

+®o)A«5(i^^;L®o)-  (4.19) 

Usually,  the  resistive  force  densities  and  energy-loss  rate  are  calculated  based  on  the  RPA 
density-density  correlation  function, 

!!<»>(?, ») 


1^2  (^5®)  I  ,  XT-rfO)/-  \  ’ 

i-v,(^)nf(^,cD) 

where  11^“^  =-{m^kgTjlT^*q)P2 ,  and 


(4.20) 


P2^°^(^,a))  =  In 


1  +  exp 


p 

g 

2m 


Y 


q  mco 


+  a 


1  -I-  exp 


p 


2m 


q  mco 
-  +  — 
2  hq 


+  a 


(4.21) 


is  the  noninteracting  electron  density-density  correlation  function,  and  v^{q)  =  e^jz  q^  is 
the  Fourier  transform  of  the  electron-electron  Coulomb  interaction  potential. 

The  scattering  of  electrons  by  impurities  such  as  ionized  donors  are  described  by  the 
screened  Coulomb  potential,  whose  Fourier  transform  is  given  by 

2 


Uiq)  = 

where  is  the  Debye  wavevector 


s(q^+ql) 


2  e^n 

q  = - 

^  e  ksTe 


fo  » 


(4.22) 


(4.23) 


with  the  temperature  and  density  dependent  factor  given  by 

f  (4.24) 

The  electron-phonon  scattering  includes  acoustic  phonons  (with  deformation  potential 
coupling  and  piezoelectric  coupling),  optical  phonons  (polar  and  nonpolar),  intervalley 
phonon  scattering,  as  well  as  interface  and  confined  phonon  scattering.  We  consider  only 
acoustic  deformation  coupling  and  nonpolar  optic  coupling  in  our  example  of  a  silicon 
device. 

In  the  case  of  deformation  potential  for  carrier  acoustic-phonon  interaction,  only  the 
longitudinal  phonon  contributes.  The  carrier-phonon  matrix  element  takes  the  form 

(4.25) 

where  d  is  the  mass  density  of  the  lattice  and  is  the  deformation  potential.  For  the 
dispersion  relation  of  the  longitudinal-acoustic  phonon  we  use  the  Debye  spectrum 

^g=^s=Vsq>  (4.26) 
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with  denoting  the  sound  speed  for  longitudinal  waves.  For  the  optic-phonon  scattering, 
we  take  the  optical  mode  frequency  as  a  constant  (the  Einstein  model): 

^gX=^o-  (4.27) 

For  the  matrix  element  of  electron-nonpolar  optical-phonon  interaction,  we  use 

(4.28) 

where  D  is  the  shift  of  the  band  edge  per  unit  relative  displacement  of  the  two  sublattices 
relative  to  the  optical  mode. 


4.3  Equivalence  to  the  Conventional  Hydrodynamic  Equations 


The  Lei-Ting  hydrodynamic  balance  equations  can  be  shown  to  be  equivalent  to  the 
hydrodynamic  equations  derived  fi-om  Boltzmann  equation  [6].  The  latters  have  the  form 


V^<|)  =  — (iV^-«), 
8 

dn  1 


dt  e 
^nw„) 

dt 


+  V.S  =  J„-E-Uw„+n{-^\, 

ot 


X  .  V)(i) = vr, 

/  KTe 


^-(^„-V)(^)  = 
ef  n 


V(|)). 


The  energy  flux  density  is 

5„=2— ^(— +-wv'), 
e  in  2 

where  Q  =  -k„V2^  is  the  heat  flux. 

The  average  single  particle  energy  is  given  by 

1  2  M 

w  =  -;rmv  H - 

2  n 


(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 


(4.34) 


and  the  energy  relaxation  rate  is  related  to  the  rate  of  change  of  the  particle  nxmiber  and 
the  energy  dissipation  rate  given  by 


,dw.  1  \u  ,dn. 


'  dt  ■ 


with 


n\n  dt' 


(-)  =-U. 


(4.35) 


(4.36) 


If  the  recombination-generation  processes  were  not  taken  into  accoimt  (C/  =  0  ),  the 
energy  relaxation  rate  reduces  to 


(4.37) 
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The  average  carrier  energy  consists  of  the  thermal  energy  and  the  kinetic  energy 

.  But  in  this  work  we  neglect  the  kinetic  energy  which  is  negligible  compared  with 
the  thermal  energy. 
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5.  The  Discretization  of  the  Semiconductor  Device  Equations 

The  hydrodynamic  balance  equations  consist  of  the  Poisson  and  the  three  balance 
equations.  These  equations  form  a  set  of  nonlinear,  coupled,  time-dependent  partial 
differential  equations.  The  solution  methods  can  be  divided  into  two  categories:  coupled 
(Newton)  or  decoupled  (Gummel)  method  [31].  The  use  of  the  Newton  method  is  limited 
by  the  requirement  of  proper  choice  of  the  initial  guess.  If  the  initial  guess  is  far  from  the 
true  solution,  convergence  is  usually  difficult  to  achieve.  Furthermore,  since  the 
equations  are  solved  simultaneously,  CPU  memory  requirement  of  the  Newton  method  is 
higher  than  that  for  a  decouple  method.  In  the  decoupled  method,  the  major  variables  are 
solved  with  currently  available  values  iteratively.  The  advantages  of  decoupled  method 
are  simplicity,  low  memory  requirements  and  convergence  for  arbitrary  initial  guesses.  In 
our  simulator,  we  use  Gummel  method  to  solve  the  hydrodynamic  balance  equations. 

It  is  well  known  that  the  numerical  solution  of  central  difference  scheme  for  the  drift- 
diffusion  equation  exhibits  oscillations  if  the  grid  is  too  coarse  or  if  the  drift  velocity  is 
too  large  relative  to  diffusion.  Although  these  oscillations  can  be  avoided  by  using 
sufficiently  fine  grids,  it  will  cost  too  much  computation  time.  This  difficulty  can  be 
circumvented  by  introducing  the  Scharfetter-Gummel  scheme  or  by  using  the  artificial 
diffiisivity.  We  chose  the  Scharfetter-Gummel  scheme  for  the  discretization  of  the  flux 
continuity  equations  in  this  work.  In  the  Scharfetter-Gummel  scheme,  a  first-order 
ordinary  differential  equation  is  integrated  by  finding  an  integration  factor.  This  approach 
results  in  an  exponentially  weighted  difference  scheme  that  prohibits  the  oscillations  in 
the  solution  without  requiring  an  excessive  number  of  grid  points. 


5.1  The  Box  Integration  Method 

Device  simulation  requires  numerical  solution  of  partial  differential  equations.  To 
solve  the  partial  differential  equations  on  a  computer,  they  must  be  discretized  on  a 
simulation  grid.  The  device  is  partitioned  into  a  finite  number  of  subdomains  as 
illustrated  in  Fig.  5.1. 
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Fig.  5.1  Illustration  of  a  typical  mesh  of  a  FET  device  for  numerical  solution  of  the 
hydrodynamic  equations. 

In  this  work,  the  spatial  discretization  is  performed  by  using  the  Box  Integration 
Method  (BIM)  [32].  The  Box  Integration  Method  is  a  generalization  of  the  finite 
difference  method.  Within  this  method,  each  node  of  the  mesh  is  surrounded  by  a 
subdomain  (box).  Each  box  is  defined  by  the  normal  bisectors  of  the  sides  emanating 
from  a  given  node,  as  shown  in  Fig.  5.2. 
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Fig.  5.2  Sample  of  Box  Integration  Method 

We  integrate  the  hydrodynamic  balance  equations  over  a  subdomain.  The  divergence 
operator  are  integrated  using  Gauss’  theorem.  So  that,  we  have  the  following  equations: 

[.p-i^dTi  =  (p-n  +  (5.1) 


f  j-i„dr^=-e\  udn. , 

Jr,  P  ”  '  Jci,  Qt  J^i 


dt 
dinw„) 


Jr,  ”  ”  '  Jn,  Qt  •’«< 


dt 


E  J„  -Uw„ -n(  "  °) 


(5.3) 

(5.4) 


f  5^  .r„dr,.  =-f  ^-dn^  +  l 

Jr,  P  ”  '  Jn,  dt  •’“< 


-  Wo 

£. j  _c/w  dn,,  (5.5) 


•/m' 


where  r,-  and  Q,  represent  the  boundary  and  the  area  of  the  /th  box,  respectively,  and  i„ 
is  the  external  vector  normal  to  r, ,  as  shown  in  Fig.  5.3. 
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The  discrete  form  of  the  RHS  of  equations  (5.1)-(5.5)  is  obtained  by  evaluating 
integrands  at  the  node  i  and  multiplying  them  by  the  box’s  area.  For  the  left  hand  sides, 
we  need  to  evaluate  the  flux  of  each  vector  through  the  boundary  of  the  box.  The  discrete 
form  of  the  LHS  can  be  expressed  as 

lFidr,='^d„F„,  (5.6) 

where  Fy  is  the  component  of  the  flux  flowing  between  node  i  enclosed  by  box  Q,.  and  a 
neighboring  node  j,  and  dy  is  the  flux  cross-section. 


5.2  Poisson  Equation 


The  electrical  displacement  vector  projected  over  the  element  side  ( i,  j  )  is: 

I'y)  .  (5-7) 

where  Ly  is  the  length  between  node  i  and  j.  Under  the  assumption  of  the  Bolt2mann 
statistics,  the  relations  of  the  electron  and  hole  concentrations  with  the  potentials  can  be 
written  as 


n  =  n.^  exp- 


.  KT,  J 


(5.8) 

(5.9) 
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where  «,>  is  the  effective  intrinsic  carrier  concentration  and\i/„and\|/p  are  the  quasi- 

Fermi  potential  of  electron  and  hole.  The  linearized  Poisson  equation  is  thus  slightly 
changed  from  its  standard  form,  and  now  reads 

2  2 

=  (5.10) 

B  e 

where  6Ay  is  defined  by 


^Dy  =  -5<|)^.) ,  (5.1 1) 

L.J 

and  5(|)  is  a  slight  deviation  of  the  electric  potential.  By  using  the  Box  Integration 
Method,  we  have  the  discretized  Poisson  equation  for  a  slight  deviation  of  the  electrical 
potential 


2^,80,. +5(1.( 


e^n 


V, 


(5.12) 


5.3  Momentum  Balance  Equation 


The  discretization  of  momentum  and  energy  balance  equations  is  carried  out  by  a 
generalization  of  the  Scharfetter-Gummel  discretization  scheme  [33,34].  The  balance 
equations  are  projected  onto  the  side  Uj  and  the  projections  are  assumed  to  be  constant 
over  the  side.  And  we  notice  that  the  first  order  differential  equation  of  the  form 


^+P(x)y  =  Q(,x), 
ax 


has  an  integral  solution  of  the  form 
y(x)  =  e  ^ 

Equation  (5.13)  can  be  rewritten  as 

>^7 


yi 


J? 


f? 


=  1. 


(5.13) 


(5.14) 


(5.15) 


JXi 

The  second  term  on  the  left-hand  side  of  (4.32),  which  is  the  convective  term,  is 
generally  small  and  is  taken  to  be  zero.  Thus,  the  electron  current  density  is  projected 
onto  the  element  side  A  leading  to  a  first  order  differential  equation 


^enksT^v^ 


dn  .1  dT^ 

+  ( 


4 


dX: 
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where  xy  is  the  coordinate  on  the  side  Ly.  In  order  to  implement  the  integration,  we  make 
an  approximation  that  electron  temperature,  electric  potential  and  total  resistive  force 
density  are  piecewise  linear  between  two  adjacent  nodes.  With  the  help  of  the 
Scharfetter-Gummel  method,  the  current  density  along  the  mesh  line  between  two 
adjacent  nodes  i  and  j  can  be  expressed  as  [21,26] 
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where  5(A)  is  the  Bernoulli  function 

A 


5(A,)^-5(-A,)^ 

^ej  ^ei 


5(A)  = 


exp(A)  - 1  ’ 


and 


T.-T. 

^ej 


k, 


■((t),-(^,)-2(7;.-r,,) 


(5.17) 


(5.18) 


(5.19) 


‘ei 

The  symbol  stands  for  the  expectation  value  of  an  arbitrary  variable  ^  over  Ly. 
Further,  we  assume  an  exponential  variation  of  the  electron  concentration  along  the  side 
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where  and  nj  are  the  electron  concentrations  at  node  i  and  j.  Thus,  the  expectation  value 
of  the  electron  concentration  is  given  by 

n  =  —  p  ndx  = - - .  (5.22) 

Ly^^i  ln(np-ln(/j;) 

The  drift  velocity  between  two  nodes  can  be  obtained  by  means  of  Jy=-enyv .  Therefore, 
the  discretized  momentum  balance  equation  is 


Similarly,  the  discretized  form  of  the  hole  current  density  is  expressed  as 
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5.4  Energy  Balance  Equation 

We  use  a  similar  procedure  for  the  energy  balance  equation.  The  electron  energy  flow 
density  is  rewritten  as 

<5-26) 

z  e 

Then  (4.3 1)  can  be  expressed  as 
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The  electron  energy  flux  density  is  projected  onto  the  element  side  I/,  leading  to  a  first 
order  differential  equation 


”  dx,  e  2  ^  ^ 


(5.28) 


Employing  the  Scharfetter-Gummel  method,  the  discretization  of  the  energy  flow  density 
can  be  written  as  [16,27,35]: 
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Finally,  the  discretized  energy  balance  equation  is 

n,  (5.32) 


dt 
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Similarly,  the  discretized  form  of  the  hole  energy  flux  density  is  expressed  as 


where 
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5.5  The  Numerical  Method  for  the  Transient  Problem 

The  time  discretization  for  the  carrier  continuity  and  the  energy  balance  equation  we 
used  here  is  the  Crank-Nicolson  implicit  scheme  [36],  which  is  second-order  accurate  in 
time  and  is  usually  stable  for  large  time  steps.  Applying  this  scheme  to  the  carrier 
continuity  and  the  energy  balance  equation,  we  have 
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Our  simulation  process  starts  with  the  potential  solution 


(5.37) 


=  --{.No  -«'^‘) .  (5.38) 

8 

While  solving  for  the  potential  at  time  level  (t+1),  the  carrier  concentration  n‘^^  at 
(t+1)  is  not  available.  So  that,  we  need  to  know  the  value  of  the  carrier  density.  If  we 
employed  the  fully  explicit  scheme  into  (5.2),  then  we  have  the  differential  form  of  the 
carrier  continuity  equation 

^?-^  =  (^V.J„-C^)'.  (5.39) 

Thus,  the  approximate  solution  of  the  carrier  concentration  at  time  level  (t+1)  is 

nf‘=M'+Ar--V-7'-Ar-C/.  (5.40) 

We  put  this  approximate  solution  into  the  Poisson  equation  to  obtain  a  modified  Poisson 
equation  [37] 

V2(|)'+‘  =  _£  (iv^  _  _  Ar .  i  V  •  r  +  AM7) ,  (5.41) 

e  e 


where  n*  and  are  the  carrier  density  and  current  density  at  time  level  (t),  respectively. 

To  avoid  the  nonconservation  of  charge  during  the  transient  state,  the  potential 
distribution  should  be  corrected  at  each  time  step  by  solving  the  Poisson  equation  [38]. 
Thus,  after  obtaining  a  solution  of  one  of  the  time-dependent  hydrodynamic  balance 
equations,  a  procedure  to  correct  the  Poisson  equation  is  necessary  to  ascertain  a  self- 
consistent  solution. 

We  will  use  the  solution  of  the  modified  Poisson  equation  to  update  the  quasi-Fermi 
potential 
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Following  that,  we  put  this  quasi-Fermi  potential  into  the  Poisson  equation  to  obtain  the 
accurate  potential  at  time  level  (t+1).  We  use  an  iterative  process  to  solve  Poisson 
equation  with  the  quasi-Fermi  potential  kept  constant  [31].  After  solving  the  Poisson 
equation,  we  have  an  accurate  potential  at  time  level  (t+1).  We  use  this  new  quasi-Fermi 
potential  and  the  new  potential  to  calculate  the  total  resistive  force  density  at  time  level 

(t+1),  ,  and  use  this  in  Eq.  (5.17)  to  get  the  current  density  at  time  level  (t+1). 

In  this  way,  we  can  solve  the  carrier  continuity  equation.  Because  the  carrier  density  has 
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already  changed,  this  caused  the  electric  potential  to  change,  too.  For  self-consistency,  we 
need  to  solve  the  Poisson  equation  again. 

The  solution  n‘'^^  is  used  to  update  the  quai-Fermi  potential  ,  which  is  given  by  [25] 


(5.43) 

Then  we  put  this  quasi-Fermi  potential  into  the  Poisson  equation  to  obtain  the  new 
potential  and  the  new  electron  density  at  time  level  (t+1).  We  use  all  these  new  solutions 


k  T 


to  calculate  the  total  resistive  force  density  and  the  energy-loss  rate  . 

With  all  these  data  at  time  level  (t+1),  we  can  solve  the  energy  balance  equation  (5.37). 
After  solving  the  energy  balance  equation,  the  electron  temperature  at  time  level  (t+1)  is 
obtained,  which  is  used  to  update  the  quasi-Fermi  potential  first.  The  updated  quasi- 
Fermi  potential  is  given  by 
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The  new  electron  temperature  will  affect  the  electron  concentration  distribution,  and  in 
turn,  the  electric  potential.  Thus,  we  need  to  solve  the  Poisson  equation  again  to  get  the 
accurate  electric  potential  and  electron  density.  Here  we  have  all  the  solutions  at  time 
level  (t+1),  which  can  be  used  as  initial  solutions  to  solve  all  the  equations  again.  The 
iterative  process  will  come  to  an  end  when  the  total  time  is  equal  to  a  predetermined 
length  of  time,  e.g.,  10  ps.  The  flow  chart  of  the  transient  simulation  procedure  is  shown 
in  the  Fig.  5.4. 

For  the  automatic  time  step  selection,  there  is  a  little  difference  between  our  method 
and  [39].  We  utilized  the  relative  difference  8*  between  the  exact  solution  and  the 
approximate  solution  of  the  carrier  density  at  time  level  (t+1)  to  select  a  new  time  step.  If 
we  take  a  time  step  and  produce  8* ,  the  new  time  step  will  be  estimated  as 


where  5  is  a  safety  factor  and  8o  is  the  given  tolerance  which  is  chosen  to  be  0.008  in 
the  following  example.  If  8*  is  larger  than  8o ,  the  new  time  step  will  decrease.  On  the 

other  hand,  when  8*  is  smaller  than  8o ,  the  time  step  size  will  increase  for  the  next 
iteration.  Special  care  is  taken  not  to  let  the  time  step  grow  too  fast.  Here,  we  restrict  the 
value  of  g' by  a  factor  1.2. 
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Fig.  5.4  Flow  chart  of  transient  solution  procedure. 
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6.  Device  Simulation  Results 


We  have  applied  our  hydrodynamic  balance  equations  to  the  simulation  of  one¬ 
dimensional  and  two-dimensional  submicron  Si  devices.  Scattering  mechanisms 
considered  are  ionized  impurity  scattering,  acoustic  phonon  scattering  via  deformation 
potential  coupling,  and  nonpolar  optical  phonon  scattering.  In  this  work,  we  included 
static  screening  only. 

6. 1  One  Dimensional  Si  n'*' -n-n**^  Diode 

Here,  a  one-dimensional  submicron  ri^  -n-ri^  diode  has  been  simulated.  It  has  some 
resemblance  to  the  n-channel  of  a  silicon  MOSFET.  Although  this  devices  is  rather 
simple,  it  is  very  useful  for  investigating  basic  transport  phenomena.  It  is  also  simple 
enough  to  do  a  Monte  Carlo  simulation,  so  that  the  latter  may  be  compared  with  our 
results.  Further,  there  are  so  many  different  device  modeling  attempts  on  such  a  device, 
so  that  we  can  readily  compare  our  results  with  those  obtained  with  different  approached. 
This  device  is  a  0.6  pm  silicon  structure  with  a  symmetric  doping  profile.  Its  inner  region 
is  0.4  pm  long  with  a  doping  level  of  Nd=2x10^^  cm'^,  while  the  emitter  and  collector  are 
each  0.1  pm  long  with  a  doping  Nd=5x10‘’  cm’^.  There  is  smooth  grading  of  the  doping 
level  at  the  jxmctions  between  the  electrodes  and  the  inner  region.  The  mesh  point  spacing 
can  be  uniform  or  nonuniform.  Here  we  used  60-node  grids  that  are  uniformly  distributed 
along  the  device  length.  The  grid  spacings  are  Ax=0.01pm.  We  chose  the  lattice 
temperature  Tq  to  be  300K. 
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Fig.  6.1  Ballistic  diode 


6.1.1  Steady  State  Simulation 

We  have  applied  the  multi-valley  Lei-Ting  balance  equation  to  a  six-valley  system 
(silicon).  The  intravalley  acoustic-phonon  scattering  and  six  phonon  models  of 
intervalley  scattering  are  considered.  Here  the  energy  band  structure  is  taken  to  be 
parabolic. 
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In  our  simulation  of  the  device,  we  have  the  left  terminal  of  the  diode  grounded  (Vl 
=0),  and  a  positive  bias  is  applied  to  the  right  terminal  (Vr=V),  so  that  the  electrons  are 
injected  fi-om  the  left  end  and  move  toward  the  right  end.  We  assume  that  the  contacts  are 
ideal  ohmic  contacts.  From  this  assumption,  it  preserves  charge  neutrality  and  thermal 
equilibrium  at  the  contacts.  We  used  the  following  boxmdary  conditions  in  this  example: 

T,{L)  =  UR)  =  T^,  (6.1) 

n(I)  =  Nd{L)  ,  n{R)  =  N^iR) ,  (6.2) 


=  (6.3) 

e  I  j  el,  ) 

where  V  is  the  applied  bias.  The  electron  velocity  and  normalized  temperature  are  shown 
in  Fig.  6.2  and  6.3.  From  our  simulation  results,  we  see  that  the  drift  velocity  exhibit 
some  overshoot  at  the  high  field  region  before  it  reaches  the  collector.  The  maximum 
drift  velocity  approaches  the  saturation  velocity  as  the  applied  bias  increases.  The 
electron  temperature  also  behaves  as  expected.  All  of  our  results  are  in  general  agreement 
with  other  available  results,  including  details  of  velocity  overshoot,  obtained  with  Monte 
Carlo  simulations  [40,41].  We  also  note  that  the  spurious  overshoot  peak,  which  usually 
exists  in  other  hydrodynamic  models  when  electric  field  drastically  decreases,  is  absent  in 
our  result. 


(cm/sec) 
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Fig.  6.2  Electron  drift  velocity  as  a  function  of  position  in  the  ballistic  diode,  under  the 
bias  voltage  of  IV,  2V,  3V,  4V  and  5V. 


Position  (mster) 


Fig.  6.3  Electron  temperature  as  a  function  of  position  in  the  ballistic  diode,  under  the 
bias  voltage  of  IV,  2V,  3 V,  4V  and  5V. 
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6.1.2  Transient  Simulation 


Here,  we  used  the  single  parabolic  hydrodynamic  balance  equation  to  perform  the 
transient  simulation. 

The  boundary  conditions  for  temperature  and  electron  concentration  on  the  emitter  and 
collector  are  the  same  as  the  steady  state  simulation.  We  applied  the  following  equations 
for  the  electrostatic  potential: 


(|)(Z;)  =  ^lnf^|  ,  ,],(/?)  = 

e  i  J  el  n,^ 


+  V(0, 


where  V(t)  is  the  applied  bias.  The  initial  conditions  were  obtained  from  the  steady  state 
solution  with  applied  voltage  V=0.  Then  at  time  t  =0,  we  applied  a  step  voltage  to  the 
anode. 


Although  the  Crank-Nicolson  scheme  is  an  implicit  time  integration  method,  it  still  has 
stability  problems  related  to  the  time  step  size[42].  In  our  simulation,  the  first  time  step  is 

l.OxlO'^^s  and  the  maximum  time  step  is  5.0xl0"'''s.  Although  we  restricted  the 
maximum  time  step  to  ensure  stability,  it  is  still  an  efficient  method  for  transient 
simulation.  For  example,  it  only  needed  221  time  steps  in  the  case  of  the  4V  bias.  From 
our  simulation  results,  we  see  that  the  drift  velocity  exhibit  some  transient  overshoot  at 
the  high  field  region  before  it  reaches  its  steady  state.  The  electron  temperature  also 
behaves  as  expected.  All  of  our  results  are  in  general  agreement  with  other  available 
results,  including  details  of  velocity  overshoot,  obtained  with  other  hydrodynamic 
simulations  [43,44]  and  Monte  Carlo  simulations  [45,46].  Velocity  overshoot  can 
improve  the  performance  of  small  device  because  the  average  carrier  velocity  can  exceed 
the  steady  state  limit. 
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Velocity  vs.  Time  (IV) 
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Fig.  6.4  (a)  Electron  drift  velocity,  (b)  Electron  temperature  as  fimctions  of  time  and 
position  along  the  device  for  a  bias  voltage  of  IV. 
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Velocity  vs.  Time  (2V) 
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Fig.  6.5  (a)  Electron  drift  velocity,  (b)  Electron  temperatiure  as  functions  of  time  and 
position  along  the  device  for  a  bias  voltage  of  2V. 
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Velocity  vs.  Time  (3V) 


Temperature  vs.  Time  (3V) 


Fig.  6.6  (a)  Electron  drift  velocity,  (b)  Electron  temperature  as  functions  of  time  and 
position  along  the  device  for  a  bias  voltage  of  3  V. 
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6.2  Two  Dimensional  MESFET 

We  have  applied  the  Lei-Ting  hydrodynamic  balance  equations  to  the  simulation  of  a 
two-dimensional  MESFET  of  the  size  0.6)mix0.2|im.  This  device  geometry  and  doping 
profile  are  depicted  in  Fig.  6.9.  The  substrate  of  the  device  is  doped  n-type  with  doping 
value  of  IxlO^^cm'^.  The  two  n'^  regions  are  of  size  O.lpmxO.OSpm.  The  doping  of  these 
regions  is  SxlO^^cm’^.  This  device  is  a  special  form  of  a  junction  field-effect  transistor. 
The  lattice  temperature  To  is  taken  as  300K. 


0.0  V  -0.8V  2.0V 


X 

n:  1.0*10*^  (cm-^)  ;  n^:  3.0*10*^  (cm'^) 


Fig.  6.9  A  two-dimensional  Silicon  MESFET  device. 
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6.2.1  Steady  State  Simulation 

We  apply  a  voltage  bias  2V  at  the  drain  and  a  negative  voltage  bias  -0.8V  at  the  gate. 
For  the  boundary  conditions,  we  assume  the  Schottky  contact  on  the  gate  and  the  ideal 
ohmic  contact  on  the  source  and  drain.  This  example  is  the  same  as  the  device  presented 
in  (Alum  et  al.  1994).  For  simplicity,  a  uniform  rectangular  grid  is  used  here. 

The  equipotential  contour  plot  is  shown  in  the  Fig.  6.10.  The  electron  concentration  is 
plotted  in  the  Fig.  6.11.  The  profile  of  the  normalized  electron  temperature  is  shown  in 
the  Fig.  6.12(a).  We  can  see  that  the  peak  temperature  is  near  the  drain.  The  longitudinal 
electron  velocity  profile  is  shown  in  the  Fig.  6.12(b).  The  spurious  velocity  overshoot 
spike  which  appears  in  the  solutions  with  the  hydrodynamic  model  (Alum  et  al.  1994)  is 
virtually  absent  in  our  results. 
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Electrostatic  Potential 
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Fig.  6. 12(a)  Electron  temperature  as  a  function  of  position. 


Fig.  6.12(b)  Electron  drift  velocity  as  a  function  of  position. 
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6.2.2  Transient  State  Simulation 

We  have  also  perfonned  transient  simulation  of  a  two-dimensional  MESFET  [28].  The 
device  is  the  same  as  before.  The  initial  conditions  were  obtained  from  the  steady  state 
solution  with  applying  a  voltage  bias  2V  at  the  drain  and  a  zero  voltage  bias  at  the  gate. 
For  the  boundary  conditions,  we  assume  the  Schottky  contact  on  the  gate  and  the  ideal 
ohmic  contact  on  the  source  and  drain.  Then  at  time  t  =0,  we  applied  a  step  negative 
voltage  bias  -0.8  V  at  the  gate. 

In  our  simulation,  the  first  time  step  is  1.0xl0''*s  and  the  maximum  time  step  is 
1.0x10"’“  s.  Although  we  restricted  the  maximum  time  step  to  ensure  stability,  it  is  still  an 
efficient  method  for  transient  simulation.  From  our  simulation  results,  we  see  that  the 
drift  velocity  exhibits  some  transient  overshoot  at  the  high  field  region  before  it  reaches 
its  steady  state.  The  electron  temperature  also  behaves  as  expected.  All  of  our  results  are 
in  general  agreement  with  other  available  results,  including  details  of  velocity  overshoot, 
obtained  with  other  hydrodynamic  simulations  and  Monte  Carlo  simulations.  The 
simulation  results  are  presented  in  the  figures  in  the  form  of  spatial  distributions  at 
various  times  during  the  entire  transient  process.  These  include  the  electron  velocity,  the 
electron  temperature  and  the  electric  field.  Also  presented  are  equipotential  contours  as 
well  as  the  terminal  currents  as  functions  of  time. 
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Fig.  6.14(a)  Electrostatic  potential  at  O.OSps 
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Fig.  7.14(b)  Electrostatic  potential  at  0.5ps 
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Fig.  7.14(c)  Electrostatic  potential  at  Ips 
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Fig.  7.14(d)  Electrostatic  potential  at  2ps 
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Fig.  7.14(e)  Electrostatic  potential  at  5.5ps 


OV  ==>  -0.8V 


0.1  0.2  0.3  0.4  0.6 

XAxis 

(micrometer) 


Fig.  7.14(f)  Electrostatic  potential  at  lOps 
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Fig.  6.15(a)  Longitudinal  velocity  at  O.OSps. 


Fig.  6.15(b)  Longitudinal  velocity  at  0.5ps. 
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Fig.  6.15(c)  Longitudinal  velocity  at  Ips. 
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Fig.  6.15(d)  Longitudinal  velocity  at  2ps. 
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Fig.  6.15(e)  Longitudinal  velocity  at  5.5ps. 


Fig.  6.15(f)  Longitudinal  velocity  at  lOps. 
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Fig.  6.16(a)  Electron  Temperature  at  O.OSps. 


Fig.  6.16(b)  Electron  Temperature  at  0.5ps. 


Fig.  6.16(d)  Electron  Temperature  at  2ps. 


Fig.  6.16(e)  Electron  Temperature  at  5.5ps. 


Fig.  6.16(f)  Electron  Temperature  at  lOps. 
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17(a)  Longitudinal  electric  field  at  O.OSps. 
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17(b)  Longitudinal  electric  field  at  0.5ps. 
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6.3  Two  Dimensional  MOSFET 

In  this  section,  we  present  our  simulation  results  of  the  behavior  of  a  deep  submicron 
n-channel  silicon  MOSFET.  This  test  device  is  the  same  as  the  device  presented  in  [48]. 
It  is  characterized  by  an  oxide  thickness  tox=7nm,  a  junction  depth  Xj=0. 1pm,  and  a 
substrate  doping  level  of  10*’cm'^.  The  length  of  the  source  and  drain  contact  is  0.15pm, 
and  the  doping  value  in  the  n^  region  is  2x10^®  cm'^.  The  effective  channel  length  is 
approximated  0.25pm.  The  device  structure  and  doping  profile  are  shown  in  Fig.  6.18.  A 
gaussian  doping  profile  is  assumed  within  the  source  and  the  drain  diffusions.  The  doping 
profile  has  been  shown  in  Fig.  6.19.  This  simulation  was  applied  on  a  21x20  nonuniform 
mesh. 

The  first  set  of  our  calculations  were  performed  on  the  silicon  MOSFET  with 
Vds=3V  and  Vgs=3V.  Fig.  6.20  shows  the  values  of  the  electrostatic  potential  for  the 
device.  The  figure  reflects  the  large  variation  of  the  electrostatic  potential  occuring  near 
the  drain  with  respect  to  the  source  and  bulk.  The  electron  concentration  and  normalized 
electron  temperature  are  shown  in  the  Fig.  6.21  and  6.22,  respectively.  The  electron 
temperature  has  a  peak  value  around  the  high  field  region  near  the  drain,  and  drops  as  the 
drain  region  is  approached.  The  electron  longitudinal  velocity  is  shown  in  Fig.  6.23.  The 
peak  velocity  is  about  1.2x10^  cm/sec,  it  is  very  close  to  the  saturation  velocity.  Fig.  6.24 
shows  the  drain  current  verse  drain  voltage  for  two  different  gate  voltages. 


Gate  3V  Oxide  thickness  =  7nm 
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Fig.  6.18  A  two-dimensional  Silicon  MOSFET  device 
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Fig.  6.20  Electrical  Potential  (Vj^g  =3V,Y^^  =3V) 


Fig.  6.21  Electron  concentration(Vjjg=3V,VQg=3V) 
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Fig.  6.22  Electron  temperature  (Vpg=3V,  V^g^SV) 
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Electron  Velocity  (Y^^=3y,y^=3y) 


.  6.23  Electron  longitudinal  velocity  (Vp5=3V,  VQg=3V) 
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7.  Discussion  and  Conclusions 

In  this  work,  a  semiconductor  device  simulator  based  on  the  Lei-Ting  hydrodynamic 
balance  equation  has  been  developed.  The  results  are  in  general  accord  with  other 
methods,  such  as  classical  hydrodynamic  models  and  Monte  Carlo  models.  However,  our 
method  treats  scattering  within  the  model  itself,  and  it  has  the  potential  of  including 
electron-electron  interaction  and  dynamic,  nonlocal  screening.  These  quantities  are 
calculated  within  the  simulation  process,  as  functions  of  the  electron  drift  velocity, 
electron  temperature,  as  well  as  the  electron  density,  without  an  outside  Monte  Carlo 
procedure.  Thus,  besides  the  usual  advantages  of  traditional  hydrodynamic  simulation 
approaches,  the  present  method  enjoys  the  added  convenience  of  self-contained  treatment 
of  scattering. 

We  have  applied  our  hydrodynamic  balance  model  to  n'^-n-n^  ballistic  diode, 
MESFET  and  MOSFET.  The  numerical  results  demonstrate  the  basic  hot  carrier  effect 
involving  the  spatial  and  transient  velocity  overshoot.  We  also  note  that  the  spurious 
overshoot  peak,  which  usually  exists  in  other  hydrodynamic  models  when  electric  field 
drastically  decreases,  does  not  show  in  our  result. 

A  generalized  Scharfetter-Gummel  discretization  scheme  with  the  Box  Integration 
method  has  been  employed  on  the  numerical  algorithm.  Also,  a  new  transient  simulation 
process  has  been  proposed  and  applied  on  the  n^-n-n^  diode  and  MESFET.  The  new 
algorithm  has  both  the  advantages  of  larger  time  steps  compared  to  conventional 
decoupled  schemes  and  lesser  memory  requirement  compared  to  coupled  scheme. 


63 


REFERENCES 

1.  G.  Baccarani  and  M.  R.  Wordeman,  “An  investigation  of  steady-state  velocity 
overshoot  in  silicon,”  Solid-State  Electron.,  vol.  28,  no.  4,  pp.  407-416,  1985. 

2.  R.  K.  Cook  and  J.  Frey,  “An  efficient  technique  for  two-dimensional  simulation  of 
velocity  overshoot  effects  in  Si  and  GaAs  devices,”  COMPEL,  vol.  1,  no.  2,  pp.  65- 
87, 1982. 

3.  M.  Lundstrom,  Fundamentals  of  Carrier  Transport,  Addison- Wesley,  1990. 

4.  K.  Hess,  Advanced  Theory  of  Semiconductor  Devices,  Prentice-Hall,  Englewood,  NJ, 
1988. 

5.  M.  Tomizawa,  K.  Tokoyama,  and  A.  Yoshi,  “Nonstationary  carrier  d3mamics  in 
quarter-micron  Si  MOSFET’s,”  IEEE  Trans.  Computer-Aided  Design,  vol.7,  no.2, 
pp.  254-262, 1988. 

6.  B.  Meinerzhagen  and  W.  L.  Engl,  “The  influence  of  the  thermal  equilibrium 
approximation  on  the  accuracy  of  classical  two-dimensional  numerical  modeling  of 
silicon  Submicrometer  MOS  transistors,”  IEEE  Trans.  Electron  Devices,  vol.  35,  no. 
5,  pp.  689-697, 1988. 

7.  C.  Jacoboni  and  L.  Reggiani,  “The  Monte  Carlo  method  for  the  solution  of  charge 
transport  in  semiconductors  with  applications  to  covalent  materials,”  Reviews  of 
Modem  Physics,  vol.  55,  pp.  645-705, 1983. 

8.  C.  Jacoboni  and  P.  Lugli,  The  Monte  Carlo  Method  for  Semiconductor  Device 
Simulation,  Springer  Verlag,  Vienna,  1989. 

9.  K.  Hess,  ed.,  Monte  Carlo  Device  Simulation:  Full  Band  and  Beyond,  Kluwer 
Academic  Publishers,  1991. 

10.  M.  V.  Fischetti,  “Monte  Carlo  simulation  of  transport  in  technologically  significant 
semiconductors  of  the  diamond  and  zinc-blende  stractures-part  I:  Homogeneous 
transport,”  IEEE  Trans.  Electron  Devices,  vol.  38,  No.  3,  pp.  634-649, 1991. 

11.  M.  V.  Fischetti  and  S.  E.  Laux,  “Monte  Carlo  simulation  of  transport  in 
technologically  significant  semiconductors  of  the  diamond  and  zinc-blende 
structures-part  H:  Submicron  MOSFET’s,”  IEEE  Trans.  Electron  Devices,  vol.  38, 
No.  3,  pp.  650-660, 1991. 

12.  R.  Stratton,  “Diffusion  of  hot  and  cold  electron  in  semiconductor  barriers,”  Phys. 
Rev.,  vol.  126,  pp.  2002-2013, 1962. 

13.  K.  Blotekjaer,  “Transport  equations  for  two-valley  semiconductors,”  IEEE  Trans. 
Electron  Devices,  vol.  17,  pp.  38-47, 1970. 

14.  G.  Baccarani  and  M.  R.  Wordeman,  “An  investigation  of  steady-state  velocity 
overshoot  in  silicon,”  Solid-State  Electron.,  vol.  28,  no.  4,  pp.  407-416,  1985. 

15.  F.  Odeh,  M.  Rudan  and  J.  White,  “Numerical  solution  of  the  hydrodynamic  model  for 
a  one-dimensional  semiconductor  device,”  COMPEL,  vol.  6,  pp.  15 1-170, 1987. 

16.  D.  Chen,  E.  C.  Kan,  U.  Ravaioli,  W.-C.  Shu,  and  R.  W.  Dutton,  “An  improved  energy 
transport  model  including  nonparabolicity  and  non-Maxwellian  distribution  effects,” 
IEEE  Electron  Device  Lett.,  vol.  13,  no.  1,  pp.  26-28, 1992. 

17.  D.  L.  Woolard,  R.  J.  Trew,  and  M.  A.  Littlejohn,  “Hydrodynamic  hot-electron 
transport  model  with  Monte-Carlo  generated  transport  parameters,”  Solid-State 
Electron.,  vol.  31,  pp.  571-574, 1988. 


64 


18.  S.-C.  Lee  and  T.-W.  Tang,  “Transport  coefficients  for  a  silicon  hydrodynamic  model 
extracted  from  inhomogeneous  Monte  Carlo  calculations,”  Solid-State  Electron.,  vol. 
35,  no.  4,  pp.561-569, 1992. 

19.  X.  L.  Lei,  J.  Cai,  and  L.  M.  Xie,  “High-field  balance  equation  for  electronic  transport 
in  weakly  nonuniform  systems,”  Phys.  Rev.  B  vol.38,  no.  2,  pp.  1529- 1532, 1987. 

20.  X.  L.  Lei  and  C.  S.  Ting,  “Green’s-function  approach  to  nonlinear  electronic  transport 
for  an  electron-impurity-phonon  system  in  a  strong  electric  field,”  Phys.  Rev.  B 
vol.32,  no.  2,pp.lll2-1132, 1985. 

21.  J.  Cai  and  H.  L.  Cui,  “Semiconductor  device  simulation  with  the  Lei-Ting  balance 
equations,”  J.  Appl.  Phys.Vol.78,  No.  11,  1995. 

22.  X.  L.  Lei,  N.  J.  M.  Horing,  and  H.  L.  Cui,  “Theory  of  negative  differential  conductive 
in  a  superlattice  miniband,”  Phys.  Rev.  Letters  vol.  66,  no.  25,  pp.3277-3280,  1991. 

23.  X.  L.  Lei,  N.  J.  M.  Horing,  H.  L.  Cui,  and  K.K.  Thomber,  “Frequency  limitations  of 
negative  differential  mobility  in  vertical  conduction  in  superlattices,”  Appl.  Phys. 
Letters  vol.  65,  no.  23,  pp.  2984-2986, 1994. 

24.  S.  Choi,  J.  G.  Ahn,  Y.  J.  Park,  H.  S.  Min,  and  C.  G.  Hwang,  “A  time  dependent 
hydrodynamic  device  simulator  SNU-2D  with  new  discretization  scheme  and 
algorithm,”  IEEE  Trans.  Computer-Aided  Design  of  Integrated  Circuits  and  Systems, 
vol.13.  No.  7,  pp.899-908,  1994. 

25.  M.  Rudan,  F.  Odeh,  and  J.  White,  “Numerical  solution  of  the  hydrodynamic  model  for 
a  one-dimensional  semiconductor  device,”  COMPEL-Int.  J.  Comp.  Math.  Electrical 
Electron  Eng.  vol.  6,  no.  3,  pp.151-170, 1987. 

26.  A.  Forghieri,  R.  Guerrieri,  P.  Ciampolini,  A.  Gnudi,  M.  Rudan  and  G.  Baccarani,  “A 
New  Discretization  Strategy  of  the  Semiconductor  Equations  Comprising  Momentum 
and  Energy  Balance,”  IEEE  Trans.  Computer-Aided  Design,  vol.  7,  pp.23 1-242,  1988. 

27.  C.  C.  Lee,  H.  L.  Cui,  J.  Cai,  E.  Lenzing,  R.  Pasture,  D.  Rhodes  and  B.  S.  Perlman, 
“Transient  device  modeling  using  the  Lei-Ting  hydrodynamic  balance  equations,”  J. 
Appl.  Phys.,  vol.  80,  no.  3,  pp.  1891-1900, 1996. 

28.  C.  C.  Lee,  H.  L.  Cui,  J.  Cai,  R.  Pasture,  D.  Woolard,  D.  Rhodes  and  B.  S.  Perlman, 
“Time  dependent  hydrodynamic  model  of  MESFET,”  Fifth  International  Workshop 
on  Computational  Electronics,  1997. 

29.  W.  Haensch  and  M.  Miura-Mattausch,  “The  hot-electron  problem  in  small 
semiconductor  device,”  J.  Appl.  Phys.  60,  pp.  650, 1986. 

30.  X.  L.  Lei,  J.  C.  Cao,  and  B.  Dong,  “Study  of  high-field  electron  transport  in 
semiconductors  using  balance  equations  for  nonparabolic  multivalley  systems,”  J. 
Appl.  Phys.,  vol.  80,  no.  3,  pp.  1504-1509,  1996. 

31.  H.  K.  Gummel,  “A  self-consistent  iterative  scheme  for  one-dimensional  steady  state 
transistor  calculations,”  IEEE  Trans.  Electron  Devices,  vol.  ED-1 1,  pp.455, 1964. 

32.  R.  S.  Varaga,  Matrix  Iterative  Analysis,  Englewood  Cliffs,  NJ:  Prentice  Hall,  1962. 

33.  D.L  Scharfetter  and  H.K.  Gummel,  “Large-Signal  Analysis  of  a  Silicon  Read  Diode 
Oscillator,”  IEEE  Trans.  Electron  Device,  vol.l6,  No.l,  pp.64-77,  1969. 

34.  A.  Forghieri,  R.  Guerrieri,  P.  Ciampolini,  A.  Gnudi,  M.  Rudan  and  G.  Baccarani,  “A 
New  Discretization  Strategy  of  the  Semiconductor  Equations  Comprising  Momentum 
and  Energy  Balance,”  IEEE  Trans.  Computer-Aided  Design,  vol.  7,  pp.231-242, 1988. 


65 


35.  A.  Leone,  A.  Gnudi  and  G.  Baccarani,  “Hydrodynamic  simulation  of  semiconductor 
device  operating  at  low  temperature,”  IEEE  Trans.  Computer-Aided  Design  of 
integrated  Circuits  and  Systems,  vol.l3,  no,  1 1,  pp.  1400- 1408, 1994. 

36.  Press,  W.  H.,  S.  A.  Teukolsky,  W.  T.  Vetterling,  B.  P.  Flannery,  Numerical  Recipes, 
Cambridge,  New  York;  Cambridge  University  Press,  1987. 

37.  S,  Yoganathan  and  S.  K.  Baneqee,  “A  new  decoupled  algorithm  for  nonstationary, 
transient  simulations  of  GaAs  MESFET’s,”  IEEE  Trans.  Electron  Dev.  Vol.39,  no.  7, 
pp.1578-1587, 1992. 

38.  M.S.  Obrecht,  M.I.  Elmasry,  “TRASIM:  Compact  and  Efficient  two-dimensional 
transient  simulator  for  arbitrary  planar  semiconductor  devices,”  IEEE 
Trans.Computer- Aided  Design  of  Integrated  Circuits  and  Systems,  vol.l4,  no.  4, 
1995. 

39.  B.  S.  Polsky  and  J.  S.  Rimshans,  “Half-implicit  difference  scheme  for  numerical 
simulation  of  transient  processes  in  semiconductor  devices,”  Solid-State  Electronics, 
vol.  29,  No.  3,  pp.  321-328, 1986. 

40.  Q.  Lin,  N.  Goldsman  and  G.-C.  Tai,  “Highly  stable  and  routinely  convergent  2- 
dimensional  hydrodynamic  device  simulation,”  Solid-State  Electronics,  vol.  37,  No. 
2,  pp.  359-371, 1994. 

41.  T.-W.  Tang,  S.  Ramaswamy,  and  J,  Nam,  “An  improved  hydrodynamic  transport 
model  for  silicon,”  IEEE  Trans.  Electron  Dev.  Vol.  40,  no,  8,  pp.  1469-1477, 1993. 

42.  P.  Wild,  “On  the  stability  of  the  time  discrtizations  for  the  semiconductor  equtions,” 
COMPEL-Int.  J.  Comp.  Math.  Electrical  Electron.  Eng. vol,  10,  no.  1,  pp.l  1-25, 1991. 

43.  E.  Fatemi,  J.  Jerome,  and  S.  Osher,  “Solution  of  the  hydrodynamic  device  model 
using  high-order  nonoscillatory  shock  capturing  algorithms,”  IEEE  Trans.  Computer- 
Aided  Design,  vol,  10,  No.  2,  pp,23 2-244, 1991. 

44.  E.  D.  Lyumkis,  B.  S.  Polsky,  A.  I.  Shur  and  P.  Visocky,  “Transient  semiconductor 
device  simulation  including  energy  balance  equation,”  COMPEL,  vol.ll.  No.  2,  pp. 
311-325, 1992. 

45.  M.  B.  Patil  and  U.  Ravaioli,  “Transient  simulation  of  semiconductor  devices  using  the 
Monte-Carlo  method,”  Solid-State  Electronics,  vol.  34,  No.  10,  pp  1029-1034,  1991. 

46.  K,  Tomizawa,  Numerical  simulation  of  submicron  semiconductor  device,  Boston: 
Artech  House,  1993. 

47.  Aluru,  N.  R,,  Kincho  H.  Law,  Peter  M.  Pinsky,  Arthur  Raefsky,  Ronald  J.  G. 
Goossens,  and  Robert  W.  Dutton,  “Space-Time  Galerkin/Least-Squares  Finite 
Element  Formulation  for  the  Hydrodynamic  Device  Equations,”  lEICE  Trans. 
Electron,  vol.  E77-C,  pp.  227-235, 1994. 

48.  W.  Lee  et  al.,  “Numerical  modeling  of  advanced  semiconductor  devices,”  IBM  J. 
Res.  Develop.,  vol.  36,  No.  2,  pp.208-230, 1992. 


